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in  a  liquid  at  low  temperatures.  Several  authors  denied 
this  possibility  and  we  show  that  such  results  are  due  to 
the  assumptions  they  made.^^ 
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part  2:  In  order  to  calculate  the  phase  diagram  of  water  we 
introduce  a  lattice  model  that  has  the  following  features. 

A  nearest  neighbor  attraction,  which  is  strongly  dependent 
on  the  relative  orientation  of  water  molecules,  due  to 
hydrogen  bonding  and  a  next  nearest  neighbor  or  three  body 
repulsion.  The  hydrogen  bonding  is  introduced  with  the  help 
of  a  set  of  weight  factors  in  accordance  with  Pauling's  ice 
rules.  The  entropy  is  calculated  according  to  the  cluster 
variation  method  for  tetrahedrons.  The  isotherms  show  a 
maximum  in  the  density  and  a  phase  separation  between  the 
vapor,  the  open  ice  state  and  a  state  which  is  dense  packed. 

Part  3:  To  see  how  orientational  dependent  forces  influence 
the  supercooling  process,  the  simplest  model  known  (Barker 
and  Fock  model)  is  analysed.  It  is  shown  under  what  conditions 
the  phase  diagram  becomes  topologically  different  from  the 
simple  Ising  phase  diagram. 

Part  4:  Preparatory  work  for  the  determination  of  the  charac¬ 
teristic  points  in  the  phase  diagram  of  an  antif erromagnet . 
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This  report  consists  of  four  parts: 

1.  Possible  Instabilities  in  the  liquid  structure  leading 
towards  the  formation  of  a  solid  (work  with  Y.  M.  Wong). 

2.  Calculation  of  a  water  model  with  the  clustervariation 
method.  Phase  diagram  and  undercooling  features. 

3.  Influence  of  orientation  on  the  undercooling  process  (work 
with  E.  Bodegom). 

4.  Transition  points  in  antiferromagnetic  systems  (work  with 


S.  Eckmecki) 
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I.  Instabilities  in  the  distribution 
function  for  a  dense  system. 

1.  Introduction. 

We  are  addressing  ourselves  to  the  question  of  whether  and 
how  the  distribution  function  of  a  dense  system  becomes  unstable 
upon  lowering  the  temperature.  As  is  evident  from  the 
observation:  a  solid  is  formed.  Just  lowering  the  temperature 
does  not  create  a  solid  however,  it  creates  a  metastable 
situation.  Does  this  metastable  situation  finally  become 
unstable  at  a  still  lower  temperature? 

Originally  Kirkwood*  pointed  out  that  this  is  the  case, 
but  in  several  subsequent  papers  it  was  pointed  out  first  by 
Kunkin  and  Frisch, 2  later  by  Nagai  and  Naitoh^  and  also  by 
Muller-Krumbhaar  and  Haus4  that  the  instability  predicted  by  the 
linearized  theory  disappeared  when  one  goes  to  higher 
approximations . 

This  note  is  to  point  out  that  this  conclusion  is  not 
necessarily  true  for  the  general  case  in  which  no  successive 
approximations  are  introduced.  Since  it  is  impossible  to 
calculate  the  general  case  completely  we  will  have  to  rely  on 
another,  but  far  more  sensible,  approximation.  This 
approximation  was  proposed  by  Kozak,  Weeks  and  Rice^  and  uses  the 
observed  fact  that  the  density  of  the  liquid  and  the  density  of 
the  solid  differ  very  little.  We  re-analyzed  the  situation  and 
and  found  that  there  is  no  contradiction  between  the  general 
liquid  theory  and  the  occurrence  of  an  instability  in  this 
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medium;  whether  such  Instability  leads  to  a  period  structure  is 
an  open  question. 

2.  The  DBGKY-formal ism. 

The  starting  equation  is  the  lowest  member  of  the  BBGKY 
hierarchy.  It  is  given  by 

g 

(See  for  instance  Balescu  )  where  ^(Rj^  is  the  singlet  density, 
u^®12^  pa*r  potential  and  g(Rj^)  the  pair  correlation 

function.  The  integral  is  over  the  volume  and  p  represents  the 
temperature.  Following  Hiroike7  we  find  from  eq.(l) 

Aa-.  C  ^  J  (2) 

V  *“■ 

where  c  ■  c(^  )  is  a  constant  in  space  depending  on  the  average 
density  ^  .  This  expression  can  be  written  in  the  Hammer stein® 
form,  if  one  puts: 

(f>  (1)  -  In  c  £(1)  (3a) 

and  - 

K(12)  >?  dr  u'(r)g(2)(r)  (3b) 

giving 

^  (1)  ♦  J  K(12)exp(^(2))d2  -  0  (4) 

Next  we  follow  the  non-linear  equation  approach  of  Week, Rice  and 


Page  4 


Kozak u  in  which  we  search  for  two  solutions  (j>  and  ^  differing 
little  in  density.  Ve  define 

(1)  -  (f>Ji  1)  (5) 
where  prefers  to  the  liquid  and  to  the  new  solution.  The 
solution  ^  satisfies  equation  (4)  with 

K(12)  -  KC 1 1  -  2 J ) , 

but  for  <£>  no  such  assumption  will  be  made.  This  is  Important 
since  this  type  of  equality  definitely  does  not  hold  for  a  solid. 
In  fact,  any  such  assumption  made  on  the  second  solution  prevents 
the  solution  of  corresponding  to  a  solid. 


differing 


Substituting  (5)  into  (4)  gives 


with 


7  (l>  ♦/ 


L(12)( 


L(12  )  =  K(12)exp^<2) 


Making  use  of  the  fact  that  the  density  of  the  new  phase  differs 
very  little  from  that  of  the  fluid  phase,  we  can  expand  the 
exponent  in  the  Kernel,  whereupon  we  obtain  the  linear  integral 
equation: 


^(D  ♦  J  L(12)*JJL(2)d2 


which  can  he  solved  by  a  Fourier  transform  in  space 


y.  <k)[l  +  L(k) ]  -  0 


For  the  existence  of  a  non-trivial  J,  { k),  one  requires  therefore 
1  ♦  L(k)  -  0  (9! 


l.e.  in  three  dimensions  one  requires 


vf 


L(R)R  sin  (kR)  dR 


(10) 


Following  Kunkin  and  Fritsch  we  will  apply  this  to  the  hard  core 
system.  This  is  not  a  very  good  system  since  it  does  not  give 
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rise  to  a  solid  in  the  ordinary  sense  of  the  word.  We  write  out 
this  example  to  show  that  one  obtains  the  same  result  in  this 
case,  but  without  the  more  stringent  assumptions  that  they  used. 


3.  Hard  sphere  system. 


For  hard  spheres,  one  has 

^  .  -i  S0'-eL'> 


(H) 


where  dQ  is  the  diameter  of  the  hard  sphere.  Inserting  this  in 


(10)  gives 


)  + 


if  J 


lc~R  =r  o 


(12) 


which  does  not  depend  on  the  temperature,  as  it  should  not.  This 
can  be  rewritten  as 

\  +.  H**?  ^  =  0  (13) 

with  z  ■  kdQ .  We  used  the  following  expression  for  the  Kernel 


L(R)  -  *(2)(do)0(do  -  R) 


4.  Discussion 

The  expression  (13)  is  exactly  the  expression  that  was 

used  by  Kirkwood,  except  that  no  inconsistent  (and  unjustified) 

linearisation  is  is  used  on  the  ^>(Rj)  and  the  g(R12),  contrary 

2  3 

to  Kunkin  and  Frisch  and  Naitoh  and  Nagai.  The  only  assumption 
is  the  step  from  (6)  to  (7),  which  is  plausible  in  the  sense  that 
it  is  in  step  with  the  observation. 
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The  occurence  of  a  second  solution,  also  called  branching 

or  bifurcation,  does  not  Imply  crystallization,  i.e.  a  periodic 

/ 

solution.  This  is  certainly  not  the  case  for  hard  spheres,  where 

g 

van  Hove  proved  ,  in  one  dimension,  that  no  periodic  solution 

exists.  For  the  non-hardsphere  case  the  question  is  left  open. 

9 

There  exists  some  work  by  Ventevogel  and  Nyboer*  on  the 
spontaneous  occurence  of  periodic  structures  at  zero 
temperatures . 

Other  approaches10  using  the  exact  solution  for  C(l,2), 
the  direct  correlation  function  or  g(l,2),  the  indirect 
correlation  functions  from  the  Percus  Yevich  equation  support  the 
absence  of  instabilities  in  three  dimension..  This  could  be  due 
to  the  Percus  YeviclV,  equation,  which  despite  its  appeal  is 
nevertheless  an  approximate  equation. 

Finally  we  remark  that  the  decoupling  schemem  used  by 

4 

Muller-Krumbhaar  and  Haus  is  not  well  understood,  that  is,  we  do 
not  know  whether  it  will  a  priori  prevent  instabll ities  to  show 
up. 
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III.  Orientational  models 


One  of  the  most  interesting  cases  in  which  orientational 

interaction  pisys  a  role  is  the  phase  diagrams  of  binary  liquids 

that  show  both  upper  and  lower  critical  temperatures.  This  phase 

diagram  was  for  the  first  time  successfully  explained  by  Barker 

and  Fock1  and  it  was  subsequently  described  in  more  detail  in  a 

2 

series  of  papers  by  Wheeler  e.a.  In  this  field  we  obtained  the 
following  results.  It  is  possible  to  map  the  Barker  and  Fock 
model  on  a  cluster  variation  type  of  description  in  a  one-to-one 
correspondence.  The  advantage  of  this  mapping  is  that  one  can 
easily  generalize  to  larger  clusters.  Barker  and  Fock  use  only 
pairclusters .  Second,  we  found  that  it  is  not  necessary  to  tie 
in  the  number  of  orientations  with  the  coordination  number  of  the 
lattice  and  one  can  generalize  the  model  accordingly. 

The  most  important  question  is,  however,  why  and  when  this 
model  does  give  the  lower  critical  point.  This  we  succeeded  in 
explaining  as  follows.  One  can  make  a  reduction,  somewhat 
similar  to  the  reductions  used  in  the  Niemei.ler  and  van  Leeuwen 
transformations  or  in  the  theories  of  decorated  lattices.  These 
transformations  lead  to  a  model  without  orientational  degree  of 
freedom  but  with  one  or  more  temperature  dependent  coupling 
constants.  By  studying  the  temperature  dependent  coupling 
constants,  mainly  graphically,  one  can  observe  under  what 
circumstances  new  critical  temperatures  will  occur. 


■>1  •"»»! 
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Usually  undercooling  can  be  obtained  very  easily  in 
substances  with  "sticky"  molecules.  Hence  we  tried  a  simple 
model  of  such  a  system.  The  model  uses  a  lattice  gas  of  two 
types  of  molecules  A  and  B  on  each  lattice  site  and  each  molecule 
has  y  contact  points.  Of  these  contact  points  y»(  are 

"sticky"  and  ^  are  not.  We  assume  that  between  the  A  and  the  B 
molecules  we  have  the  following  energy  expression:  if  either  one 
or  both  contact  points  of  a  pair  of  molecules  are  sticky  they 
attract  each  other  and  if  both  are  of  the  non-sticky  kind  they 

\ 

repel.  We  choose  this  example  since  Barker  and  Fock  have  shown 
that  this  case  leads  to  a  lower  critical  point.  However,  it  is 
easy  to  show  that  this  is  one  example  out  of  a  class  of  systems. 

The  entropy  for  such  a  system  in  the  notation  of  Kikuchi 
is  given  by  (we  start  out  with  the  pair  aporoximation) : 

s;  o  TTU.LM*'' 

^  “  TT  tyoj  L)  !  3,;i  L  J 

where  L  is  the  number  of  lattice  points,  xA  and  xB  the  site 
variables  and  y^j  the  sixteen  pair  variables  of  the  problem.  The 
coordination  number  is  z  and  the  weight  factors  g^j  are  given  by: 
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the  resulting  free  energy  per  site  is  given  by 


f  J  ~  -  S‘J  y«j‘  _l<1 

~  *  Vo-  y*a- ^ a“')V 


where  z  is  the  number  of  nearest  neighbors. 


The  constaints  on  the  pair  variables  are 


*  “  f  ^  i 

b-  a  y>j  i  »«•-** 


with  the  appropriate  Lagrange  multipliers  behind  the  semicol* 


With  these  terms  inserted,  we  minimize  the  extended  free 


energy 


f  * kT  1 S  s  ( **  -  j  y*  ^ 


leads  to 


y.«  ”  e  *  >  y;4-  =  cxf.{(^f7jVj^ .  & 

since  we  assumed  no  interaction  between  like  molecules. 


-P?M 


The  chemical  potential  Ka  is  found  from  n>*  a 

which  Rives 

Pa  =  -UT^-O^iaX^  +  lUr^-V^C) 

= -Ut  |o-o  xA  + 1 

and  similar  for  If  the  solution  has  two  roots  (PA)j  and 

(^A)2  we  deal  with  coexisting:  phases. 

To  solve  the  system  we  simplify  the  notation.  Define 

=  ^  Mr6"1**0 

We  find  the  followinp  system  of  equations  for  the  variables 

Xj  ,  •  •  .x^ 

X;  2  Xj  ^  *• 

where  xj  “  x2  "  xa  and  x3  "  x4  "  XB* 
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IV.  Antiferromagnetic  transition  points. 


We  address  ourselves  to  the  question  of  the  transition 

points  of  a  mixed  antiferro  and  ferromagnetic  system;  a 

so-called  metamagnet.  The  global  phase  diagram  of  a  metamagnet 

consists  of  a  Neel  temperature  (a  transition  to  the 

antiferromagnetic  state  at  zero  field)  and  a  line  of  second  order 

transitions,  a  possible  tri-critical  point,  a  line  of  first  order 

transitions  for  zero  staggered  field  and  lines  of  first  order 

transitions  in  non-zero  staggered  fields.  These  options  were 

described  by  Kincaid  and  Cohen1  who  showed  how  the  various  phases 

can  be  obtained  from  a  Landau-Ginsburg  expression  in  the  various 

order  parameters.  If  we  ignore  the  non-analytic  behavior  near 

the  critical  points,  this  free  energy  is  an  analytic  expression 

and  they  computed  the  various  coefficients  using  the  mean  field 

theory.  Under  this  assumption  they  find  that  there  are  two 

regions  for  the  coupling  constant  ratio  J/J  =  5  .  One  region 

F  AF 

2.  ^  2.  *  leads  to  a  tri-critical  point,  the  other  to  a  so-called 
bi-critical  end  (BC)  point.  We  found  2)  that  the  next  higher 
approximation,  the  Bethe  approximation,  leads  to  a  set  of  phase 
diagrams  with  BC  points  only.  Despite  the  fact  that  this 
approximation  gives  generally  much  better  results  for  the 
(absolute)  critical  temperature,  there  are  reasons  to  believe 
that  the  Bethe  approximation  does  not  take  enough  detail  into 
account.  (By  absolute  critical  temperature  we  mean  to  refer  to 
those  experiments  where  both  the  temperature  and  the  coupling 
constants  are  known  and  not,  as  is  unfortunately  always  almost 
the  case,  where  the  critical  temperature  is  measured  without  any 
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further  Information  about  the  coupling  constants.)  Returning  to 

the  fact  that  the  Bethe  approximation  may  not  always  Rive  the 

details  of  thephase  diagram;  a  most  striking  example  is  the  work 

3 

on  the  AB  alloys  where  only  the  tetrahedron  approximation 
started  to  give  a  complete  phase  diagram.  Consequently  we  want 
to  study  the  metamagnet  using  the  cluster  variation  method. 

The  cluster  variation  method  requires  more  information 
about  the  lattice  than  the  Bethe  approximation.  In  the  latter  it 
suffices  to  specify  the  number  of  nearest  and  next-nearest 
neighbors.  For  larger  clusters  we  have  to  give  the  detailed 
information  not  only  of  the  lattice,  but  also  of  the 
substructures  we  expect  to  obtain  in  the  antiferromagnetic 
ordering.  After  some  discussion  we  decided  to  study  the 
following  cases:  1)  the  bcc  subdivided  in  two  sc,  2)  the  sc 
subdivided  in  two  fee,  and  3)  the  fee  in  subdivided  two 
substructures  as  mentioned  in  Li.^  Finally  we  evaluated  the  2 
dim.  square  into  two  2  dim.  square  lattices  in  order  to  have  an 
test  case. 

At  this  moment  we  have  determined  the  free  energy 
expressions  for  these  four  cases.  In  order  to  calculate  the  free 
energy  in  the  clustervariation  method  one  needs  the 
clusterrelations  between  the  various  order  parameters  and  their 
weight  factors  in  order  to  construct  the  corresponding  entropy 
expressions.  This  work  is  being  done  by  Mr.  Servet  Eckmecki . 
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ABSTRACT 

In  order  to  calculate  the  phase  diagram  of  water  we 
introduce  a  lattice  model  that  has  the  following  features.  A 
nearest  neighbor  attraction,  which  is  strongly  dependent  on  the 
relative  orientation  of  water  molecules,  due  to  hydrogen  bonding 
and  a  next  nearest  neighbor  or  three  body  repulsion.  The 
hydrogen  bonding  is  introduced  with  the  help  of  a  set  of  weight 
factors  in  accordance  with  Pauling's  ice  rules.  The  entropy  is 
calculated  according  to  the  cluster  variation  method  for 
tetrahedrons.  The  isotherms  show  a  maximum  in  the  density  and  a 
phase  separation  between  the  vapor,  the  open  ice  state  and  a 
state  which  is  dense  packed. 


Keywords:  water,  equation  of  state,  cluster  variation  method, 

entropy,  ice,  maximum  density  in  water,  hydrogen  bond. 
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INTRODUCTION 


1  2 

In  two  of  his  recent  papers,  Pell  ’  developed  a  lattice 
model  in  calculating  the  phase  diagram  of  water.  He  used  n  bcc 
lattice  and  placed  oxygen  atoms  and  vacancies  on  the  lattice 
points;  hydrogen  atoms  in  water  molecules  lie  on  bcxiy -diagonal 
bonds  connecting  nearest-neighbor  lattice  points. 

The  bcc  lattice  is  made  of  two  interpenetrating 

diamond-type  sublattices,  and  thus  we  can  define  the  ordered  and 

disordered  phases  based  on  these  two  sublnttices.  In  this  way  of 

treatment,  liquid  water  is  in  a  disordered  state  that  lies 

between  two  ordered  states.  One  ordered  state  is  the  dense  bcc 

state,  associated  with  ice  VI,  and  the  other  state  has  only  half 

of  the  available  sites  occupied  (i.e.  one  sublattice  is  almost 

fully  occupied  and  the  other  sublattice  almost  empty).  This  open 

structure  we  associate  with  ice  I.  Which  of  the  states  is  formed 

depends  of  course  on  the  pressure  and  the  temperature.  Although 

for  the  full  treatment  of  these  two  phase  transitions  one  needs  a 

two  sublattice  model,  in  this  introductory  Part  I  we  will  explain 

the  key  points  of  the  treatment  without  using  the  sublattices  in 

order  not  to  "omplicate  mathematics.  The  sublattice  treatment  is 

3 

presented  in  the  accompanying  Part  II. 


In  order  to  make 

the 

lattice 

model 

as  water-like 

as 

possible,  we  build  in 

two 

features , 

as  was 

done  by  Bel 1 . * 

The 

first 

feature  is  the  presence 

of  the  hydrogen 

bond  between 

two 

water 

molecules.  The 

hydrogen  bond 

has 

the  result  that 

two 

neighboring  molecules  have  an  interaction  potential  that  depends 
on  the  relative  orientation:  if  mutually  aligned  in  the  proper 
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way,  there  is  a  very  strong  binding  force;  if  not  mutually 
aligned,  the  force  between  the  two  water  molecules  is  weak.  If 
the  water  molecule,  which  has  a  V-shape,  is  placed  in  the 
body-center  of  a  cube,  the  two  hydrogen  "arms"  can  be  laid  in  12 
different  ways  along  the  eight  body  diagonal-halves  that  radiate 
from  the  body-center.  It  is  tacitly  assumed  that  the  water 
molecule  always  prefers  solely  these  orientations,  even  in  the 
absence  of  neighbors.  If  we  pick  a  pair  of  water  molecules,  the 
number  of  orientations  will  be  144,  and  it  is  easy  to  show  that 
18  of  these  lead  to  hydrogen  bonding.  The  combinatorial 

4 

calculation  is  in  fact  identical  with  the  one  made  by  Pauling  to 
obtain  the  entropy  of  ice  at  zero  degrees. 

The  second  feature  of  the  Bell  model  is  the  use  of  a 
repulsive  three-body  force.  Bell  reasons  that  this  force 
discourages  too  close  a  packing  and  may  be  the  main  reason  why  a 
negative  expansion  coefficient  in  water  is  found.  Whether  an 
actual  three-body  force  really  exists  is  an  open  question, 
despite  the  fact  that  such  forces  are  often  proclaimed  in  the 

5 

literature,  because  there  is  really  no  direct  evidence  for  its 
presence  in  nature.  All  conclusions  were  based  on  the  necessity 
to  fit  the  data.  The  lattice  models  work  .  h  the  same  kind  of 
handicap:  except  for  the  hydrogen  bond  energy,  we  do  not  have  an 
independent  and  consistent  estimate  of  the  other  coupling 
parameters  in  our  model ;  hence  we  have  to  work  backwards  from 
the  known  isotherms.  In  this  paper  we  restrict  ourselves  to  work 
with  one  set  of  parameters,  taken  from  a  melange  of  information 
available . 
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Figure  1  shows  the  lattice  structure  we  use  in  our 
formulation.  The  positions  A  form  a  diamond  lattice  and  the 
positions  B,  which  form  by  themselves  also  a  diamond  lattice, 
will  complete  this  to  a  hcc  lattice.  In  the  treatment  in  this 

paper  we  choose  a  tetrahedron  as  indicated  in  Fig.  1  as  the 

basic  cluster.  In  this  tetrahedron  two  of  its  corners  are  on  the 

A  sublattice  and  two  of  its  remaining  corners  are  on  the  B 

sublattice.  Of  the  six  edges,  four  are  nearest  neighbors  and  two 
are  next-nearest  neighbors. 

The  method  used  to  obtain  the  expression  for  the  free 

T  B 

energy  was  developed  by  one  of  us  ’  *  and  was  successfully 

applied  to  a  large  variety  of  models.  In  Section  2  we  will  give 

a  description,  but  the  main  idea  is  the  following.  If  one  writes 

down  the  entropy  associated  with  a  cluster  of,  say  four, 

particles  using  the  Boltzmann  expression  one  has  made  an 

overcount.  This  overcount  has  to  be  corrected  by  subtracting 

partial  entropy  contribution  due  to  clusters  of  smaller  size. 

The  entropy  of  these  subclusters  again  contains  an  overcount, 

which  has  to  be  corrected  by  considering  still  smaller  clusters, 
q 

and  so  on."  The  result  differs  substantially  from  the  expressions 
used  by  Bell  and  his  co-workers.  In  his  one-sublattice  case*  he 
introduces  the  entropy  of  the  tetrahedrons  and  of  the  point 
variables,  but  omits  all  intermediate  clusters.  On  the  other 
hand  in  the  two-sublattice  paper  Bell  and  Salt^  use  the 
point-variable  entropy  only,  which  reduces  the  model  to  the 
simplest  mean  field  calculation. 
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In  order  to  facilitate  the  computations  we  introduce  a 
cluster  relation  CR-matrix,  given  in  Appendix  I.  Using  this 
CR-matrix  we  can  write  down  the  expression  for  the  entropy  and 
hence  the  grand-potential  function  ^  .  Minimization  cf  the 

latter  leads  then  to  the  most  probable,  i.e.  equilibrium, 
distribution  of  the  largest  clusters.  This  variation  of  the 
with  respect  to  the  order  parameters,  ten  in  our  case,  leads 
to  a  set  of  self-consistent  equations.  Kikuchi  has  pointed 
out10,11  that  this  set  of  equations  can  be  solved  in  a  natural 
way.  Natural  meaning  that  the  corresponding  free  energy  is 
always  decreasing  during  the  iterative  process.  This  method  is 
used  here,  too,  and  the  solutions,  where  obtianed,  are 

practically  at  the  limit  of  the  accuracy  of  the  computer. 

2.  Description  of  the  model . 

In  Figure  2  we  give  the  set  of  ten  clusters  introduced  by 
Rell  as  well  as  all  the  resulting  subclusters.  The  open  circles 
refer  to  the  absence  of  a  particle,  the  closed  circles  to  the 
presence  of  an  oxygen  atom.  Hydrogen  atoms  attached  to  the 
lxygen  atom  are  not  shown  except  for  the  11-bond  cases.  One 
connecting  line  represents  a  nearest  neighbor  bond,  two 
connecting  lines  represent  a  next-nearest  neighbor  bond.  The 
nearest  neighbor  bond  is  either  "blank"  or  has  an  "H"  (actually 
two  crosslines)  to  indicate  whether  that  particular  pair  of 
occupied  sites  is  or  is  not  hydrogen-bonded.  Note  that  there  are 
two  types  of  pair  bonds:  y  and  z.  The  z-bonds  are  next-nearest 
neighbor  bonds  and  do  not  depend  on  the  relative  orientation  of 
the  two  water  molecules  at  the  pair,  while  y-pair  bond  refers  to 
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noarcRt  neighbors  i\ has  two  options:  with  and  without  an 

11 -bond  ,  resulting  in  the  Paul ing-rat io  of  weight  factors:  one  to 

seven.  In  the  table  wo  distinguish  between  the  orientational 

weight,  factors  g  on  the  right,  and  the  topologicnl  weight 

factors  g^  on  the  left.  The  g^  is  the  number  of  different 

configurations  when  the  molecules  and  h-bonds  are  moved  around  in 

the  tetrahedron.  The  g  is  the  number  of  wavs  the  hydrogen  atoms 

o 

can  be  placed  next  to  the  oxygen  molecules.  The  two  weight 

factors,  the  topological  and  the  orientational,  enter  the 
calculation  in  a  different  way.  In  what  follows,  anything  that 

is  simply  called  the  weight  factor  is  the  product  of  the 
two:  g  ^  =  (R^g0)j  for  the  i  *  configuration. 

In  order  to  derive  the  cluster  relations  we  observe  that 
each  cluster  can  be  augmented  by  one  more  site  to  form  the  next 

larger  cluster.  Take  for  instance  te  cluster  called  u j .  If  this 

cluster  is  completed  into  a  tetrahedron,  the  added  site  can  be 
either  occupied  or  not  occupied.  This  leads  to: 

tlj  *  Wj  +  12w^  (2.1) 

since  occupation  can  take  place'  in  12  different  orientational 
ways  of  the  hydrogen  atoms. The  parameters  x  can  be  completed  into 
a  linear  combination  of  z’s  and  also  into  a  linear  combination  of 
y’s.  In  such  cases  both  relations  should  hold.  This  procedure 
is  not  unambi guous .  There  is  a  bifurcation  in  y0  which  can  be 
written  as  two  different  linear  combinations  of  u*s: 

y 2  "  u;»  +  12u4  "  Ug  +  ipu,.  +  gun  (2.2) 

These  subclusters  are  defined  in  Figure  2.  It  is  clear  that 
during  the  construction  of  these  relations,  one  should  use  the 
orientational  weight  factors.  The  normalization  relations  of  the 
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w's,  tho  u's,  the  z's,  the  y's,  and  tho  x's,  iiro  tho  total  weight, 
factors.  Tho  result  of  tho  calculation  Is  Riven  In  Appendix  I, 
in  tho  form  of  a  matrix  whore  all  triple,  pair  and  single  cluster 
parameters  are  written  as  linear  combinations  of  the  tetrahedral 
cluster  parameters,  w^  (i  «  1,...,10).  In  Table  I  the  y^ 
expression,  for  example,  satisfies  both  the  reght-hand  u  sums 
(2.2).  We  now  introduce  the  our  model. 


3.  Pe terminat  ion  o_f  the  entropy  . 


The  entropy  for  the  tetrahedral  model  on 
cubic  lattice  can  be  written  symbolically  ns  : 


a  bodycentered 


S  = 


(3.1) 


where  N  is  the  total  number  of  lattice  points  in  the  system. 

This  entropy  expression  is  known  to  be  the  most  efficient  one 

7  8  10  11 

when  the  tetrahedron  is  used  as  the  variahle  ’  ’  ’  for  bcc . 

The  powers  in  (3.1)  are  derived  from  the  procedure  of  correcting 

over-counting  as  mentioned  in  section  1.  The  curly  bracket 

7  8  10  11 

notation  is  the  standard  one  used  in  the  CVM  '  *  ’  and  is,  for 

example, 


N  SIV  (Nx^!*! 

Using  Stirlings  ap?>roximatlon , 


(3.2) 

we  can  write  (3.1)  explicitly  as 
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S/kK-  -  fc  'i  fix  | 


(3.3) 


.00 


where  of  (x)  is  short  for  (x)rr  xlnx  -  x,  nnd  whore  the 

variables  u,  z,  y,  and  x  are  to  be  expressed  In  the  variable  w. 

We  would  like  to  point  out  that  in  Doll's  calculation1  the  term 

in  w  and  the  term  in  x  wore  included  but  that  the  terms  in  u,  z, 

12  13 

and  y  wore  omitted.  Our  expression  (3.3)  is  known  '  *  to  Rive 
more  reliable  results  than  the  one  used  by  Roll.  Wo  first  make 
the  formulation  symmetric;  for  instance  wo  write: 

3£(*ij)  -  3/2(«((ziJ)  +*C(*kl))  (3.4) 


and  then  compute  the  derivative  which  can  be  written  in  general 
as : 


whereiB^^  are  the  elements  of  the  decomposition  matrix  (P-matrix) 
which  is  to  he  explained  in  (3.6)  below. 

Let  us  make  a  short  remark  about  the  new  labels:  the 

4 

subscripts  "i.lkl"  represent  2  -  16  options.  In  a  non-hydrogen 

bonded  model  they  can  be  replaced  by  five  options  with 
appropriate  weight  factors  g^.  In  Fig.  2  they  have  to  be 
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replaced  by  ten  options  to  distinguish  between  the  possible 
presences  and  absences  of  the  hydrogen  bonds.  We  think  the 
figure  speaks  for  itself. 

The  elements  of  the  decomposition  matrix  can  be  written 
down  immediately  by  considering  the  number  of  ways  a  tetrahedral 
cluster  can  be  broken  up  into  its  constituent  elements;  that  is, 
into  triangles,  double  bonded  pairs,  single  bonded  pairs,  and 
sites.  We  found  a  simple  relation  between  the  elements  of  the 
CR-matrix  and  the  elements  of  the  D-matrix: 


(3.6) 


where  fj^  equals  four,  except  for  the  parameters  z,  where  f^  =  2. 
This  relation  also  has  the  practical  value  that  it  simplifies  the 
data  input  in  the  computation. 

The  breakdown  of  the  tetrahedrons  into  subclusters  is  not 
entirely  unambiguous .  The  point  variables  xg  represent  an 
occupied  site.  The  molecule  on  such  an  occupied  site  can  be 
oriented  in  two  different  sets  of  6  directions  and  this  will 
determine  how  the  diagram  has  to  be  completed  into  a  y-pair. 
Each  completion  can  lead  to  a  different  linear  combination  of 
pair-clusters  using  y^  or  y4  and  there  is  no  a  priori  guarantee 
that  these  two  equations  give  the  same  result  for  the  variable 
Xg.  The  correct  way  to  assure  this  is  to  introduce  a  separate 
Lagrange  multiplier,  which  leads  to  a  minor  iteration  in  the 
natural  iteration  method.11  An  alternate  way  is  simply  to  take 


L . . . - 
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the  average,  in  this  case  one  to  seven.  We  computed  both  cases 
and  found  not  much  difference  in  the  final  result  for  the  order 
parameters  of  the  10  tetrahedrons,  at  least  not  in  the  region  of 
temperatures  and  chemical  potential  used.  The  issue  was  not 
further  persued  because  in  the  extended  model  the  two 
orientations  of  the  molecule  on  a  given  site  were  explicitly 
introduced  in  the  model  and  consequently  this  bifurcation  does 
not  occur . 

To  construct  G/N,  the  grand  potential  per  lattice  site,  we 
have  to  add  two  more  terms  to  the  expression  for  TS.  One 
represents  the  energy  of  the  various  occupations  of  the  cluster 
and  one  term  contains  the  chemical  potential;  in  total: 

£  2  (1;  -  F/a.v  )  3<-'uV 

(3.7) 

Where  ai  is  the  number  of  occupied  sites  in  the  cluster  w^ ,  i.e. 

the  number  of  black  circles  in  w^^  of  Fig.  2.  The  expression  for 

the  energies  are  written  as  linear  combinations  of  the 

three  basic  parameters:  the  pair  energy  in  the  absence  of  a 

/r*  ‘ 

hydrogen  bond,  ■£&  £  the  additional  pair  energy  in  the  presence  of 

H 

a  hydrogen  bond  and  ^  the  next  nearest  neighbor  repulsion. 

The  numerical  factors  stem  from  the  fact  that  there  are  6 
tetrahedrons  per  lattice  site,  containing  a  total  of  24  sites. 

Finally  we  add  a  Lagrange  multiplier  term,  of  the  form: 

-  Z 

t's(  V 


(3.8) 
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In  order  to  maintain  the  normalization  condition. 

Combining  these  relations,  the  grand  potential  G  is 
written  as 

-  6^ t  (v - 3,^6«0 


When  this  is  minimized  with  respect  to  w^  ,  we  obtain,  for 
example , 


+  LLtii,  -  lL(n 


(3. 10) 

-  p>Y"\=  c, 


This  expression  has  an  easily  understandable  physical  meaning. 
When  we  look  at  w^  in  Fig.  2,  we  see  that  it  contains  2u2’s, 
2ug's,  2z2's,  yj,  2y2's,  y^,  2xj's  and  2x2's.  These  subclusters 
appear  in  the  equation.  The  coefficients  on  logarithm  terms 
originate  in  Eq.(3.9).  All  of  the  derivatives 
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(*  “  1,  ..  .  .10),  which  arc  derived  from  Fq.(3.9),  can 


be  interpreted  in  the  same  way. 


When  these  derivatives  vanish,  we  can  see  the  meaning  of 
by  forming 


§  »  $  -  f.  *w.-  *  BX  -  p  £ 

1  i= i  i  r 


(3.11) 


Hence  the  Lagrange  multiplier  X  is  equal  to  G/N  when 


iteration  has  converged. 


From  the  grand  potential  we  find  immediately  the  pressure 
p,  since  it  is  given  by 


"V  G  U  V  L  «■’ 


(3. 12) 


where  a  is  the  length  of  the  edge  of  the  bcc  cube.  In  order  to 
determine  the  equation  of  state  we  need  to  know  the  density  . 
This  quantity  can  bo  either  expressed  in  terms  of  the  a^s  that 
are  used  in  the  chemical  potential  term  as  : 


Z.  'U>;  =  <a> 

iti.  ^ 


(3. 13a) 


or  directly  by  using: 


izXi 


(3. 1 ;3b) 


the  relation  between  these  two  expressions  can  he  obtained  by 
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writing  x0  In  torms  of  the  order  parameters  as  given  by  the 
column  of  the  CR-matrix.  We  find 


lo 

2  C. 


(3.14) 


Next  we  can  compute  the  two  additional  observable  quantities: 
the  number  of  nearest  neighbor  pairs,  18  (7y_  +  y  )  and  the 
number  of  next  nearest  neighbor  pairs,  144  z„.  These  quantities 
have  been  reported  in  X-ray  experiments  and  have  the  fascinating 
property  that  they  remain  rather  "ice-like"  in  the  low 
temperature  region  of  the  liquid.*^ 


Finally  we  report  a  similar,  but  not  directly  observable, 
quantity:  the  relative  number  of  hydrogen  bonded  nearest 
neighbor  pairs: 


y-i 

7  + y* 


(3. 15) 


Although  this  model  cannot  give  the  nngular  correlation  between 
two  neighboring  molecules,  R^  is  in  a  certain  way  a  measure  of 
this  angular  correlation.  The  latter  can  be  experimentally 
extracted  from  the  results  of  polarized  light  scattering 
experiments. 


4.  Choice  of  F.nergy  Parameters . 


Pago  lf> 


The  model  wo  used  contains  four  adjustable  parameters. 
One  more,  the  lattice  constant,  is  needed  in  case  the  grand 
partition  function  is  converted  into  a  pressure.  These  four 
parameters  are: 

1.  the  nearest  neighbor  attractive  interaction  in  the 
absence  of  the  hydrogen  bond , 

2.  the  additional  attractive  energy  introduced  when  the 
relative  orientation  between  a  pair  of  molecules  leads 
to  a  hydrogen  bond, 

3.  a  repulsive  energy  that  discourages  the  simultaneous 
occupation  of  next  nearest  neighbor  sites.  This  can 
be  introduced  in  the  form  of  a  repulsive  three  body 


force  or 

through 

two 

body 

' next-nearest 

neighbor ) 

repulsive 

force . 

The 

latter 

seems  to 

us  more 

realistic . 

It  has  been  suggested  that  there  may  be  a  three  body  force 
detectable  in  the  vapor  phase  of  water  and  other  gases.  Direct 
evidence  for  such  a  force  is  hard  to  come  by.  ‘Very  often  the 
three  body  force  was  inserted  in  the  calculations  of  the  virial 
coefficients  to  make  up  for  discrepancies  between  the  computation 
and  the  experiment.  The  presence  of  open  ice  suggests  that  in  a 
certain  range  of  temperatures  and  pressures,  the  occupation  of 
next-nearest  neighbors  is  discouraged.  (We  can  justify  this  in  a 
schematic  way  by  proving  that  the  number  of  relative  orientatons 
with  repulsive  energy  exceeds  the  number  of  pair  orientations 
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with  attractive  energy  for  next-nearest  neighbor  sites;  see 
Appendix  II.)  This  can  be  accomplished  in  the  model  by 
introducing  either  a  next-nearest  neighbor  repulsive  force  or  by 
a  repulsive  three  body  force. 

As  to  the  most  realistic  values  for  the  interaction 

constants,  even  if  the  interaction  between  water  molecules  were 

precisely  known,  ’  ’  it  is  hard  to  assess  how  potential 

energy  functions  would  have  to  be  translated  into  the 

parameters  of  such  a  schematic  model.  Roughly  speaking  the 

interaction  energy  £^is  about  1.5  kcal/mol  (corresponds  to 

li/k  =  700K).  To  estimate  the  H-bond  potential  we  use  the 

16 

probability  of  broken  H-bonds  as  introduced  by  Luck  and  Ditter  . 
Their  result  for  this  probability  P(T)  can  be  represented  by 

TCT)  •*-*)?  (  Hy«5/r) 

which  leads  to  H/k  =  -500  for  the  H-bond  potential.  Although 
these  values  should  not  be  taken  too  seriously,  it  is  interesting 
to  notice  that  they  lead  to  a  critical  temperature  (and  critical 
pressure)  that  lies  in  the  neighborhood  of  the  observed  values. 


The  next-nearest  neighbor  repulsion  is  very 
assess.  Hence  we  are  only  guided  by  the  need  that 
lie  above  (2/3) $  in  order  to  obtain  open  ice  at 


val ues . 


hv\ 


difficult  to 
its  value  must 
low  pressure 


5.  Results 
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Tho  computations  wore  executed  for  a  number  of  potentials. 
We  report  here  only  one  or  two  cast's,  since  it  became  clear  that 
this  simple  model  always  lacks  details  at  high  densities  since 
the  average  over  two  orientations  were  taken.  The  main  search 
was  to  see  whether  we  could  reproduce  the  density  maximum  that 
was  found  in  the  molecular  field  model.  The  potential  that 
showed  the  desired  features  most  clearly  is  given  by  =  1000, 

=  2500,  and  *2.  =  500  K.  See  Figure  3.  The  isotherms  are 

n 

given  for  five  different  temperatures:  T  =  300,  325,350,  375  and 
400  K.  The  points  were  obtained  by  lowerinp  the  absolute  value 
of  JJ>  ,  the  chemical  potential.  Jl)  is  nepative.  In  each  new 
point  the  values  of  w's  of  the  preceding  point  are  used  as 
initial  conditions  for  the  iteration.  Each  isotherm  starts  at  a 
high  density-high  pressure  point  and  then  comes  down  in  pressure. 
At  a  certain  value  of  the  chemical  potential,  and  corresponding 
pressure,  the  number  of  iterations  increases  and  the  calculation 
becomes  susceptible  to  a  transition  to  a  different  value  of  the 
density.  In  this  particular  calculation  this  happens  twice, 
except  for  the  high  temperature  isotherm  which  is  smoothly  going 
from  high  density  to  low  density.  It  can  be  seen  from  the 
chemical  potential  versus  ^  at  low  temperatures,  that  for  the 
given  potential  there  are  two  different  slopes,  one  between  0  and 
1/2  and  a  different  one  between  1/2  and  1,  consequently  there  are 
two  transitions:  one  from  the  low  density  or  gas  phase  to  the 
one-half  density  or  open  ice  phase  and  one  from  that  phase  to  the 
high  density  or  liquid  phase. 
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Tho  con tor  part  of  tho  plot  shows  that  the  Isotherms  are 
crossing  each  other.  This  Implies  that  there  is  a  density 

maximum,  since  the  density  is  the  same  at  two  different 

temperatures.  This  is  plotted  separately  in  Figure  4.  The 
dotted  line  indicates  the  phase  transition. 
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Cluster  relation  (OR)  matrix. 


This  matrix  represents  the  value  of  each  subeluster  as 
expressed  in  a  linear  combination  of  the  principal  clusters.  The 
10  principal  clusters  are  the  possible  occupations  of  the  four 
sites  of  the  tetrahedron  with  and  without  hydrogen  bonds. 
Although  this  information  can  be  deduced  from  the  text,  we  repeat 
it  here  so  that  we  can  use  it  in  future  extensions  of  the  model. 
The  matrix  C  is  related  to  the  decomposition  matrix  D  in  a  simple 
manner,  see  (3.G).  We  used  this  relation  to  check  the  matrix 
elements  of  C.  The  results  are  given  in  Table  I. 
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Appendix  II 


Next  nearest  neighbor  repulsion. 


Since  the  next-nearest  neighbor  repulsion  seems  to  be  an 

essential  part  of  the  model,  we  give  the  following  suggestion 

Suppose  we  deal  with  3molecules,  as  shown  in  Fig.  Al;  each 

molecule  has  a  dipole  moment,  and  takes  12  orientations  dictated 

— ^ 

by  the  model.  Then  the  dipole  moment  will  be  oriented 

(twice)  in  the  ♦  x,  *y ,  and  +z  directions,  and  the  dipolar 
coupling  energy  between  1  and  2  and  between  2  and  3  is  given  by: 

-  e  (pj^y  -  px  p!  -  p*  p* )  j 

using  the  dipolar  tensor.  If  we  consider  the  six  possible 

— > 

directions  of  |J(  ,  we  find  that  the  lowest  energy  is  obtained  for 

24  different  arrangements.  If  we  now  consider  the  relative 
orientations  of  dipole  1  with  respect  to  dipole  3  we  find  that  16 
have  zero  energy,  4  have  energies  that  mutually  cancel  out  and  4 
are  repulsive.  These  4  are  antiparallel  along  the  z  axis. 

We  may  repeat  this  combinatorial  excercise  assuring  the 
presence  of  an  11-bond  between  either  1  and  2  or  2  and  3  and  we 
find  again  that  the  relative  orientations  between  2  and  3  in 
which  the  dipolar  interaction  is  repulsive  are  favored.  Since 
this  reasoning  depended  on  the  presence  of  an  intermediary 
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moloculo,  this  model  scorns  to  favor  the  concept  of  an  effective 
three  body  force. 
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List  of  Figures 


1.  Insertion  of  the  tetrahedron  in  the  bee  lattice  and  its 
subdivision  in  two  diamond  lattices. 

2.  Tetrahedrons  with  and  without  hydroqen  bonding,  and  the 
resulting  subclusters  required  to  obtain  the  Kikuchi  entropy 
expression. 


3.  Isotherms  for  T  =  300,  323, 

330, 

373  and 
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potential  piven  by  2  =  - 

1300, 

i?  =  2300 

and  ST  =  300  K. 

The 

H 

crossinp  of  the  isotherms  menus  a 

density 

maximum . 

4.  Direct  computation  of  the  density  maximum  for  the  same 
potential  . 

5.  (Al).  Positioning  of  three  molecules:  1-2  and  2-3  are 
nearest  neighbors  and  1-3  are  next-nearest,  neighbors. 
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